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I. INTRODUCTION 



The propagation of pulses in granular materials has been a subject of vigorous recent interest. Seminal work 
on this subject was carried out by Nesterenko ^ two decades ago in his analysis of the propagation of nonlinear 
compression pulses along a line of particles. This early work firmly established the nonlinear flavor of the problem: 
Nesterenko showed that under appropriate assumptions, among them the slow spatial variation of the displacements 
of the particles, the equations of motion for granular particles could in most cases be approximated by a continuous 
nonlinearpartial differential equation that admits a soliton solution (later shown to actually be a solitary ^^^^ 
' n\ r. agating perturbation in the chain. The recent revival of interest in the subject [3, 's', 4,'?,'^ 




has been tri gger ed in part by a concern with important technological applications 
such as the design of efficient shock absorbers |14(. the detection of buried objects [^0,0,13, and the fragmentation 



of granular chains The revival has involved advances in the modeling, simulatioUjand solution of the equations 
associated with important features of granular materials such as their discreteness 0, S ILCLi IM i dimensionality , 
disorder 0, 0, Q| , and loading provided by gravitational forces 0, 0, 0, 0, 0, 0| ■ The preponderance of the work 
has been numerical, but important analytic advances have also been made, as well as experimental verifications of 
some of the theoretical predictions 6, 16]. 

The standard generic model potential between monodisperse elastic granules that repel upon overlap according to 
the Hertz law is given by [l^, 123 

F(4,fe+i) =a\5\l,^,, S<0, 

(1) 

V{5k,k+i) -0, S>0. 

Here 

Sk,k+i = 2R- [{zk+i + Vk+i) - {zk + Vk)] , (2) 



"^^"^ D{Y^a)^l(^±—^], (3) 



5L'(y, cr) V 2 y ' 2 V y 

Y and a denote Young's modulus and Poisson's ratio, Zk denotes the initial equilibrium position of granule k in the 
chain, and yk is the displacement of granule k from this equilibrium position. The geometric parameter R is the 
radius if the particles are elastic spheres. More generally, R is determined by the principal radii of curvature of the 
surfaces at the point of contact [2(|. The exponent n is 5/2 for spheres, it is 2 for cylinders, and in general depends 
on geometry. 

The displacement of the fc-th granule (fc = 1,2, . . . ,_L) in the chain from its equilibrium position in a frictional 
medium is governed by the equation of motion 

^~dT^ = ~ "(y*^ ~ yk+iT^^d{vk - vk+i) + a{yk-i - vkT^^oivk-i - vk)- (4) 



Here 9{y) is the Heavyside function, 9{y) = 1 for y > 0, 0{y) ~ for y < 0, and 9{Q) — 1/2. It ensures that the 
particles interact only when in contact. Note that for open boundary conditions the second term on the right hand 
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side of this equation is absent for the last granule and the third term is absent for the first, while for periodic boundary 
conditions y^+L = Uk- In terms of the rescaled variables and parameters 

/ 9\l/n -, / o\l/n ~ / 9\l/n 

\ a J vq \ a J mvo \ a J 

Eq. (0J can be rewritten as 

Xk = -JXk - {Xk - Xk+l)'^~^0{xk - Xk+l) + (.Tfc-l - Xk)"~^0{xk-1 - Xk), (6) 

where a dot denotes a derivative with respect to t. 

Initially the granules are placed along a line so that they just touch their neighbors in their equilibrium positions 
(no precompression), and all but one granule, granule i, are at rest. The velocity of granule i is t^o (the impulse). In 
the rescaled variables the initial conditions become Xk{0) = ife(O) — Q,yk ^ i, Xi{0) = 0, and Xi{0) = 1. In our work 
we set i — 1. 

In the absence of the frictional contribution —"fdyk/dr or —jxk, when n > 2 an initial impulse settles into a 
pulse that becomes increasingly narrow with increasing n, and propagates at a velocity that is essentially constant 
and determined by n and by the amplitude of the pulse. For n — 2 the pulse spreads in time and travels at a 
constant velocity independent of pulse amplitude. In the latter case there is considerable backscatteringthat leads to 
backward motion of all the granules behind the pulse, whereas the backscattering is minimal for n > 2 The pulse 
is a completely conservative solitary wave in the limit n — > cx). 

Our interest lies in ascertaining the effects of the frictional contributions on these results. Although friction is on 
occasion mentioned in theoretical and numerical work, it is usually mentioned as an element that is neglected or 
omitted. However, its presence and importance in experiments is recognized as inevitable 0. Note that frictional 
effects may arise not only from the forces between the granular chain and the surrounding medium but also from the 
conversion of translational motion to other degrees of freedom (e.g. rotation) |16| |. In an earlier paper we analyzed 
in detail the effects of frictional forces on the dynamics of two granules, specifically the way in which forward and 
backward motion of the granules is affected by these forces Herein we extend that work to a chain of granules. 

Since our approach will in general follow established methods, it is useful to lay out at the outset an overview 
of the principal approaches that have been implemented in the study of pulse dynamics in frictionless chains, and 
the various features that determine these dynamics. Theoretical studies of pulse dynamics in frictionless chains have 
mainly relied on three approaches: 

(1) Numerical solution of the equations of motion [IItL Igl lTollllllllll^ : 

(2) Continuum approximations to the equations of motion followed by exact or approximate solutions of these 
approximate equations P, 0, ; and 

(3) Phenomenology about properties of pairwise (or at times three-body) collisions together with the assumption 
that pulses are sufficiently narrow to be principally determined by these properties 0, |23| . 

From these studies it is clear that there three features determine the dynamics in these chains: 

(1) The power n in the potential; 

(2) The absence of a restoring force; and 

(3) The discreteness of the system. 

While each of these leads to essential aspects of the dynamics, the literature is not always entirely clear on which 
feature is determinant in a particular behavior, nor is it always clear which of these features is being approximated 
or even ignored. One example is the equivocal connection between the dynamics of granular systems and systems 
with power law interaction potentials that include a restoring force. For instance, there is an extensive literature on 
the Fermi-Pasta-Ulam (FPU) chain with purely nonlinear interactions of the form appearing in Eq. ^ with n > 2 
but without the Heavyside thcta functions . Highly localized pulses propagate in these systems (l^ |2^ , but their 
relation to the localized solutions in granular systems is by no means clear. This ambiguity is exacerbated by the fact 
that discrete systems are frequently approximated by continuum equations. While the pulse-like (soliton or solitary 
wave) solutions that emerge from these approximations are assumed (and even shown numerically) to describe certain 
aspects of granular system, the continuum approximations never explicitly respect the absence of a restoring force, 
i.e., they are more clearly justifiable for FPU-like systems. And so with these solutions in hand, it is not clear what 
consequences of the absence of a restoring force and of the discreteness of the system have been lost. In some sense 
opposite to the continuum approximations are two-body or three-body phenomenologies. These provide considerable 
qualitative insight into some consequences of discreteness, but are quantitatively not sufficiently accurate because 
the true dynamics extend over more than two or three granules. In this paper, even while focusing on the effects of 
friction, we make an attempt to provide some clarity on these issues. 

Numerical simulations are of course extremely helpful and ultimately an accurate test of theory. Note that the 
rescaled system Eq. (0 has no free parameters in the absence of friction and only a single parameter when there is 
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friction, and therefore numerical characterization is particularly straightforward in this particular system. However, 
this universality is lost with any of a number of modifications that might be interesting and have been considered 
recently such as, for example, tapering of the masses in the chain 0,0,0] and mass 0,llll or frictional disorder. 
Consequently, it is desirable to have a strong analytic backdrop. 

The paper is organized as follows. Section ^] deals with a chain of cylindrical granules. First we review existing 
results for frictionless granules, discussing them in the context of the issues and uncertainties noted earlier. We then 
extend the theory to granules subject to friction, and present numerical simulations in support of these results. In 
Sec, mil we follow the same presentation sequence for granules with n > 2, with special attention to spherical granules 
in our numerical simulations. A summary of results and of future directions is presented in Sec. II VI 

II. CYLINDERS (n = 2) 

Cylindrical granules have provided a theoretical testbed for dynamics in granular chains because the exponent 
n = 2 leads to analytic manageability not available for other potentials. Although sometimes called the "harmonic" 
case, it should be remembered that the Hertz potential is quite different because there is no restoring force, that 
is, the cylindrical Hertz potential is half of a harmonic potential. The derivative of the force law therefore changes 
discontinuously between extension and compression. This leads to considerable differences in the chain dynamics. 
Nevertheless, some aspects of the cylindrical chain dynamics can be inferred from those of a harmonic chain, and this 
can be exploited to great advantage in the analysis. 

If one were to ignore the absence of a restoring force and implement the simplest continuum approximation to 
Eq. in the absence of friction, the result would be a simple wave equation with diffusive coupling whose solutions 
do not represent the observed behavior of the n = 2 chain. In reality, an initial impulse in the chain described by 
Eq. © in the absence of friction moves as a spreading pulse. Although the pulse spreads and sheds some energy, the 
waveform can nevertheless be clearly identified as a pulse Q . Its maximum kmax (t) travels forward with a constant 
unit velocity and a displacement amplitude that increases as t^^^. The pulse spreads in time as t^/^, more slowly 
than it would in a system with diffusive coupling. It is interesting to understand which of these features are due to 
the discrete nature of the chain, and which are due to the absence of a restoring force. Further, once these questions 
have been answered, it is interesting to explore what happens to these features in the presence of friction. 

Our presentation of the chain of cylindrical granules consists of three parts. First we review the results of Hinch 
and Saint-Jean (HSJ) 9] for a frictionless chain, and we recast the problem in a way that clarifies some of the 
approximations made in that work and pinpoints the sources of differences between the n = 2 chain and a purely 
diffusive coupling. We supplement this review with our analytic results that reproduce some of their purely numerical 
ones. Then we modify this analysis to include the effects of weak hydrodynamic friction on the granules. Finally, 
we complement this analytic (and necessarily approximate) treatment with a comparison with numerical simulation 
results for the frictional chain of cylindrical granules. 

A. Frictionless Granules - Theory 

The HSJ theory is based on the following approximations implemented consecutively and independently. 

1. The solution is assumed to be described by a traveling pulse of constant form which propagates at constant unit 
speed and has an amplitude b and width A that vary slowly with time: 

x(fc,t) = 6W/(^). (7) 

2. The pulse is assumed to retain almost all of its initial energy. HSJ show that assumption 1. leads to equipartition 
of this energy between potential and kinetic forms, as should be the case for a harmonic potential. 

These two assumptions are sufficient to lead to the conclusion that A oc 6^. 

3. A continuum approximation is implemented that takes into account some discreteness effects. In the lowest 
order diffusive coupling approximation one would set the difference Xk+i + Xk-i — 2xk ~ d^x(k)/dk^ . Retention of 
the next term in a Taylor series expansion of Xk±i about Xk leads to the continuum approximation that incorporates 
some of the effects of discreteness: 

d^x{k,t) _ d^x{k,t) 1 d^x{k,t) 

~~d^ dk^ ^12 9fc4 ■ ^' 
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Note that this expansion includes a restoring force. A transformation to a moving frame with unit propagation velocity, 
i.e., a change of variables from k and t to v = k — t and i, transforms this equation to 

d'^x d'^x 1 d'^x 

(9) 



4 ■ 



5*2 dtdiy 12 9j/ 

Regardless of the form of b{t) or of the function /, this is sufficient to establish that the solution Eq. Q is consistent 
with this equation only if A ~ t^^^, and also that the next term in the Taylor series expansion is unimportant for this 
result. The width of the pulse is therefore governed by the first manifestation of discreteness. One assumes that the 
absence of a restoring force does not affect this result because the granules within the pulse are in fact overlapping most 
of the time and hence subject mainly to the repulsive portion of the potential. Note that the continuum approximation 
together with the conservation of energy are then sufficient to arrive at the conclusion that b{t) ~ t^^^, i.e., that the 
pulse amplitude actually grows with time. 

4. It is not yet clear that (|7|) is actually compatible with (jSJ until one determines the function /. Substitution of 
(0) into and retention of leading terms in t leads to the equation for f(£,), 

/"" - 8e/" - 4/' = 0. (10) 

This equation has four solutions. The one consistent with the requirement that f{£^) decays for large (i.e., ahead of 
the wave) and consistent with the assumption of energy conservation by the pulse is (26| 

/•oo 

fiO^N I M'{2'/'y)dy (11) 

where Ai{z) is an Airy function and 



N : 



Er, 



/^7^z4(2l/3y)dy 



1/2 

= (12) 



Here Ep is the pulse energy and 2^^^£^q is the first zero of Ai{z). 

These features describe the traveling_displacement pulse of increasing width and amplitude quite accurately, as 
shown by the numerical simulations in 9] . We note that the total energy of the system with the initial unit velocity 
impulse at one granule is 1/2. The numerical results of HSJ lead to an asymptotic pulse energy of Ep = 0.48, that is, 
96.2% of the energy resides in the pulse. Below we calculate the pulse energy analytically. 

In addition to the forward traveling pulse, HSJ also investigate the momentum of the particles ejected backward as 
the pulse goes by. That particles must be ejected is a consequence of the conservation of momentum: the traveling 
pulse of constant energy and increasing width and displacement amplitude carries increasing forward momentum, 
which must be balanced by the backward momentum of the ejected particles. One then arrives at the next item on 
the list of assumptions: 

5. The absence of a restoring force is explicitly recognized in the calculation of the momentum of the ejected 
particles, which simply keep traveling backward with the momentum they acquire at the moment of separation from 
the pulse. Equating the rate of change of the momentum of the forward pulse at time t to that of the particle ejected 
at that time leads to the conclusion that the backward momentum of the n"* particle 

the latter equality arising from the unit speed of pulse propagation. The numerical simulations of HSJ lead to the 
value c = 0.158, a value that we obtain analytically. 

To obtain analytic results for the pulse energy Ep and the constant of proportionality c, we note that the forward 
momentum of the propagating pulse is 



m I nm = mt) (13) 



p= E W^^) 

x{k,t)>0 x{k,t)>0 ^ ' 



?0 



where 



/'OC 

N^N Ai^ {2^ ^ ^ ^d^ ^1. 379 ^/E^. (14) 
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The rate of change of the forward momentum then is P — Nb{t). Since the pulse velocity is unity, this is the 
momentum transferred to the last particle as it is ejected: 



The backscattered energy then is 



Vb^-P = -Nb(t) = = _:?^„-5/6. (15) 

6 6 



(16) 




and for the total energy we obtain 

E^Ep + Eb = 1.056Ep, (17) 

from which it immediately follows that asymptotically the energy in the pulse is 

Ep = 0.947^;. (18) 

Thus, 94.7% of the energy resides in the pulse (to be compared to the HSJ value of 96.2% obtained from simulations) 
and the remainder is backscattered. With the initial condition E — 1/2 used in all simulations we thus have for the 
constants defined earlier N = 2.431 and N — 0.949. The constant c in HSJ is c = A^/6 = 0.158, exactly as they 
obtained from numerical simulations. 

This essentially completes the solution. The summary description is then that an initial velocity impulse propagates 
forward at unit speed, with a width that grows as t^^'^ and a displacement amplitude that grows as 

^1/6 This pulse 

carries almost all of the initial energy and its momentum increases. This increase in momentum is compensated by 
granules that are ejected backward as the pulse passes. The speed of the ejected granules decreases, the n*'' granule 
being ejected with a speed proportional to 

Further insights can be gained by viewing this problem a bit differently. Suppose that we do not implement a 
continuum approximation at all, but instead focus on the full discrete equation 

Xk ^ Xk+i + Xk-i - 2xk. (19) 

Like the continuum equation, this of course includes restoring forces. This equation can be solved exactly [23 |. 

Xk{t)^ f J2k-2{2t')dt' , (20) 

where Jn{z) is the Bessel function of the first kind '2^. The solution is a moving spreading pulse along with oscillatory 
displacements that are left in its wake (precisely because there are restoring forces) . The peak and width of the pulse 
can be obtained from the properties of the Bessel functions, in particular from knowledge of their first two zeroes and 
the first maximum [23|. The maximum of the pulse occurs at kmax{t) = t + 0{t^/^), so the pulse velocity is unity. Its 
width increases as i^/'^. In these features the solution is similar to the one assumed in Eq. Q. However, the amplitude 
of the displacement pulse does not grow but is instead constant in time. (If we write the Bessel function solution 
in the form Q but with b{t) = const independent of t, we find from substitution into Eq. ||HJ) that / satisfies the 
equation / — 8z/ —8f— 0.) The pulse energy according to this description is not constant but instead decreases 
as t^^^^. The energy that is lost goes into the oscillatory displacements in the wake of the harmonic pulse caused by 
the restoring forces. The additional and appropriate assumption of energy conservation in HSJ adds this lost energy 
back into the pulse without affecting its spreading rate or velocity. This indicates that (in the n = 2 problem) energy 
conservation is an appropriate additional assumption in lieu of the absence of a restoring force. 

In summary: the fact that the solution of the n — 2 Hertz problem is a spreading pulse of unit velocity and of width 
that increases as i^/'^ is purely a consequence of the discreteness of the system and not dependent on the presence or 
absence of restoring forces; it is a feature of a discrete harmonic system. The (near) conservation of energy in the 
pulse is not a feature of a harmonic system and must be included as an additional assumption. Conservation of energy 
must be implemented not only to describe all the features of the traveling pulse but also to describe the backward 
momentum of the particles ejected as the pulse moves along. Extensive numerical results quantitatively supporting 
the features just described can be found in 0. 
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B. Granules With Friction - Theory 



Here we generalize the previous theories to a chain in which the granules experience friction [cf. Eq. ((HJ]. We 
generalize the HSJ theory to this case and also use the solution of the damped harmonic chain to complement these 
results. 

1. The solution is assumed to be a traveling pulse of constant shape which propagates at unit speed, has a width 
A that varies slowly in time, and an amplitude that aside from an exponential decay due to the friction also varies 
slowly with time: 



x{k, t) 



- p-7t/2 



mf 



k - t 



(21) 



We retain only the 7-independent portions of the slowly varying functions A and b. This is thus a "lowest order" 
ansatz. Indeed, we will show later that corrections to this lowest order ansatz seem to be of 0(7^). The pulse energy 
then decays exponentially as e~^*' . Explicitly, the kinetic energy is obtained as follows. Retaining only leading orders 
in time, the pulse velocity that follows from Eq. H21|) is 



x(k, t) 

Therefore, the kinetic energy to 0{'^t) is 



m 



0(7t). 



(22) 



b{t) 



f 



2 X{t) 



?0 



•b'it) 



2 A(t) 



where Ep is the undamped pulse energy. On the other hand, the potential energy may be written as 

2 



U{t) = \Y,[xk+i-Xkf 



1 



2I: 

k 



dx{k, t) 



dk 



which also leads to 



U{t) 



A(t) 



to 0{'^t). The total pulse energy to this order then is 

Ep{t) -- 



x(t) ■ 



(23) 



(24) 



(25) 



(26) 



2. We assume that almost all of the energy resides in the decaying pulse. 

These results and assumptions are thus again sufficient to conclude that A oc 6^. It then follows from Eq. (|26|l that 
the only effect of friction to this order is thus the overall exponential decay of the energy. 

3. As before, a continuum approximation is implemented that takes into account some discreteness effects and of 
course the frictional contribution (while still including the nonexistent restoring force, as before). For this purpose 
it is convenient to implement the change of variables z(k,t) = e'''*/^x(fc, i). Retention of the next term in a Taylor 
expansion of Zk±i about Zk beyond the purely diffusive approximation leads to 



d^z{k,t) d^z{k,t) 1 d^z{k,t) Y 



dt^ 9fc2 

The associated moving frame equation [v = k — t) is 

d'^z d'^z 



dt^ dtdv 



12 9fc4 

1 d'^z 
129^^ 



z{k,t). 



With the ansatz 



z{k,t)^h{t)f 



(27) 



(28) 



(29) 
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associated with (|21|) this is again sufficient to establish that the pulse widens as A '-^ t^^^. Corrections begin to set in 
when the last term in Eq. H28|) becomes important, i.e., for times t > 7"^/^. For small 7 this is a time much longer 
than the entire lifetime of the pulse, which is of 0(7^^) due to the overall exponential decay. As before, we are again 
led to the conclusion that b{t) ~ t^^^. 

4. The envelope function is again the solution of Eq. IjlOl) with corrections of 0(7^). 

5. As before, as the pulse travels forward granules are ejected backwards. In the frictionless problem these granules 
continue moving backward with the same velocity forever after ejection because they are simply freely moving granules. 
In the problem with friction these granules progressively slow down because they, too, are damped, but it is of interest 
to calculate their momentum at the moment of ejection from the pulse. Indeed, we now show that in the presence 
of damping the backward ejection momentum is greater than in the undamped chain, a result consistent with that 
found in a two-granule system |2ll | . 

The forward momentum of the propagating pulse to 0{"/t) is 

x{k,t)>0 x(k,t)>Q ^ ' 



= e~^"%{t) / /'(Od^ = Ne-^'/^it), (30) 

Ho 

where N is given in Eq. (|14|l . The rate rate of change of the forward momentum is 



V 2b(t)J 



(31) 



The last particle in the pulse loses momentum 6 through dissipation as the pulse moves across it. This momentum 
loss is the difference in the momentum of this particle when it is in the middle of the pulse (i.e., when it is maximally 
compressed, at which point its momentum is zero) and when it is at the end of the pulse and about to be ejected (at 
which point its compression is zero): 

S = -jAx = --fxmax - -7iVe-T*/26(t). (32) 
Since the pulse velocity is still unity, the momentum transferred to the last particle as it is ejected is 



Vb = -P 
N 



V 2b{t) bit) J 



— e 
6 



2b{t) b{t) 



To facilitate comparison of this result with that of the frictionless chain we expand the exponential and retain terms 
of 0(7i): 

vb^-^t-'^''(l + \^ty (34) 

This confirms that the backward ejection momentum is indeed greater than in the undamped chain. 

The forward moving pulse of the n = 2 Hertz problem with small friction is essentially identical in shape to that 
of the frictionless problem except for an overall exponential decay. The pulse travels at unit speed and its width 
increases as t^/^ ^ these two features again being a consequence of discreteness and essentially unaffected by friction. 
The assumption that the pulse energy decreases only because of the friction is an additional assumption. This feature 
must again be implemented separately to describe the amplitude of the traveling pulse correctly, and also to describe 
the backward momentum of the granules ejected as the pulse moves along. 

An analysis of the full discrete equation (with restoring forces) starts from the linear equation 

Xk^Xk+i+Xk-i-2xk-'-iXk- (35) 
The change of variables Zk{t) = e^*/^Xfe(i) immediately leads to 

72 

Zfe = Zk+i + Zk-i - 2zk + —Zk (36) 
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FIG. 1: Total energy decay as a function of time. The energy decays exponentially as e (lines). The symbols are the 
numerical simulations: circles for 7 — 0.001, squares for 7 = 0.005, diamonds for 7 = 0.007. 



and consequently to the solution 



Xk{t) 



J2k-2i2t')dt' + Oi-f^) 



(37) 



The solution is a moving spreading decaying pulse along with decaying oscillatory displacements. As before, the 
peak of the pulse occurs at kmax{t) = t + 0{t^^'^), its width increases as t^^"^, and the amplitude of the pulse decays 
exponentially. 

In summary: the fact that the solution of the n = 2 Hertz problem with friction is an exponentially decaying 
spreading pulse of unit velocity and of width that increases as t^^^ is purely a consequence of the discreteness of the 
system and not dependent on the presence or absence of restoring forces; it is a feature of a discrete harmonic system. 
The assumption that the pulse retains almost all of the system energy (which decays exponentially) is again not a 
feature of a harmonic system and must be included as an additional assumption. This conservation feature must be 
implemented not only to describe correctly all the features of the traveling pulse but also to describe the backward 
momentum of the particles ejected as the pulse moves along. This backward momentum is greater in the chain with 
friction than in the frictionless system. 



C. Numerical Simulations 



In this section we check the accuracy of our analytic results for the n = 2 frictional chain through numerical 
simulations. First, we note that the total energy of the system decreases exponentially, as seen in Fig. ^ for various 
values of the friction coefficient. For times t < 107~^, aside from the overall exponential decay of the energy the 
system behavior is indeed the same as that of the frictionless system in that it broadens as t^^^ and propagates at 
unit velocity. We find that the pulse ceases to exist at a time of order t ^ 107""'^. More specifically, we find that at 
a time ^ 8.67"^ the pulse energy decreases abruptly (more rapidly than exponential), and at a time ^ 15.77^^ the 
backscattered energy becomes greater than the pulse energy. 

We must check our prediction for the self-similar impulse wave /(^) propagating along the chain of particles. We 
follow HSJ and plot the scaled velocity as a function of the scaled position for different times in Fig. [3 From Eq. H22|l 
the appropriate scaled velocity to leading order in the friction is /'(C) = xXe^*'/'^ /b oc if^/^e^*/^, and this is the 
ordinate in the panels. The abcissa is ^ = (fc — i — 0.25)/A oc {k — t — 0.25)i~^/^. The shift 0.25, which also occurs in 
the frictionless case, is carefully explained and derived in HSJ and comes about because the scaling function solution 
has velocities in the propagating impulse wave of 0{t'^^^^) (here modified by the exponential friction factor) while the 
rebound velocity of particles is of 0{t^"^^^) (again modified by the exponential factor). This indicates that there must 
be a correction of C'(t~^/'^) to the leading-order term, and this correction appears as a shift in the scaling. The solid 
curve is /' as obtained from Eq. while the points are the results of our numerical simulations. The agreement is 
very good provided the damping is small, but serious deviations begin to set in with increasing damping, as seen in 
panel (d). 
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FIG. 2: Chain of cylindrical granules: scaled velocity pulse as a function of scaled position for different friction coefficients. 

We have not calculated further corrections analytically, but a numerical test is possible. One might test a correction 
of 0(7) [cf. Eq. 123)], e.g. of the form b{t) = t^/^ [1 + -fC{i^,t)], or of 0(7^) [cf. Eq. The coefficient C in the 

correction would in general be a function of v. We have tested various corrections with the (admittedly unjustified) 
assumption that C does not depend on v (the numerical effort to do otherwise seems unwarranted), and have found 
that the leading correction seems to be of 0(7^), specifically b{t) = t^^^{l+27j^t). This result is analytically appealing 
but is purely a numerical outcome. The rescaled velocity results are shown in Fig. |3| 

Finally, we have predicted a somewhat unexpected feature of the backscattered granules, namely, that at the 
moment of ejection their velocity is greater than that of the corresponding ejected granules in the absence of friction 
[cf. Eq. H34|) ]. This is seen in Fig.^ In the first panel we show the velocity of the last ejected granule as a function of 
time for different values of the friction coefficient. The velocity is indeed higher with increasing friction, up to a time 
beyond which more complex behavior not captured by our lowest order theory sets in. In the second panel we have 
scaled the velocities by the factor (1 + ^jt) which, according to our theory, should collapse the curves. At very early 
times, while the collapse occurs with the scaling predicted by our theory, the resulting curve is not yet the theoretical 
one because the pulse is not yet clearly defined. However, after this early period we see both the collapse as well 
as agreement with the theoretical line up to times at which higher order effects set in, indicating that our theory 
captures the correct behavior up to those times. 



III. SPHERES (n = 5/2) 



The archetypal system for the study of impulse dynamics associated with a Hertz potential of exponent n > 2 
consists of spherical granules. Such a nonlinear potential gives rise to a traveling pulse of essentially constant speed 
that seems to retain its shape for very long times, i.e. it travels with essentially constant amplitude and constant 
width. We say "essentially" because this is the outcome of approximate theories J_, 9] and of numerical simulations 
but is not a rigorous result (except for n — > 00). The pulse arises from the balance of the dispersive forces that tend 
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-2-10123 -2-1012 
(x-t-0.25)t"'"^ (x-t-0.25)t"'"^ 

(c) 7 = 0.010 (d) 7 = 0.020 



FIG. 3: Chain of cylindrical granules: velocity pulse corrected for higher order frictional effects as a function of scaled poisition 
for different friction coefficients. 




log(t) 



(a) Unsealed (b) Sealed 



FIG. 4: Backward velocity of the last ejected granule as a function of time. The curves are for 7 = (thick line), 0.001 
(dashed-dotted line), 0.005 (dashed line), 0.007 (dotted line), and 0.01 (thin line). The line with the star symbols in the second 
panel has slope —5/6. The slope of the frictionless curve is exactly —5/6, as predicted by HSJ. 
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to spread the excitation, and the nonlinear and discrete nature of the system that focuses it. As we saw in the last 
section, discreteness is not sufficient to stop the widening of the pulse, but only to slow it down. Nonlinearity is 
necessary to obtain a pulse of constant width, at the very least approximately. The result is a narrow pulse, involving 
no more than a handful of particles at any time. The pulse is narrower with increasing n, but it is already sharply 
defined in the classic case n — 5/2. This sort of behavior is also known from the classic FPU problem where the 
potential typically considered is quartic, n — A (or a combination of n = 4 and lower order contributions), and includes 
restoring forces |24| . The absence of dispersive restoring forces in a Hertz system causes the localized excitation to be 
even more stable than in the associated FPU system. 

In this section our analysis again consists of three parts. We begin by briefly reviewing the theory for a chain of 
frictionless spherical granules as first presented by Nesterenko ^] . Then we modify this analysis to include the effects 
of friction. In this case there is no exactly solvable discrete counterpart even if one includes restoring forces. Finally, 
we present numerical results to support our findings. 



A. Frictionless Granules - Theory 



The theory first introduced by Nesterenko 1] and later augmented and complemented by others is based on the 
same approximation implemented in the case of cylindrical granules, i.e., a Taylor expansion of Xk±i about Xk and 
retention of a term beyond the first (which here is no longer purely diffusive coupling). This continuum equation, 
assumed to hold within the compression pulse, again does not explicitly exclude restoring forces and is therefore 
equivalent to a continuum approximation with some contributions due to discreteness for the FPU problem. The 
starting equation is thus 



Xk = -{xk - a;fc+i)" ^ + {xk-i 



Xk) 



and the resulting continuum approximation is 



d'^x 



d_ 

dk 



dx 
"dk 



24 



dx 
dk 



^ d^x 


1 d^ 


Y dxy-'' 


dk^ 


~ 2Adk^ 


\ dk) 



(38) 



(39) 



This equation reproduces Eq. ^ when n = 2. 

As before, one implements a change of variables to a moving frame, with ^ — k ~ cot. Here cq is a speed to be 
determined. The big difference between the equation one obtains with n > 2 and that obtained earlier for n = 2 
[Eq. ©J is that it only involves the variable the variable t no longer appears explicitly: 



d_ 

d^ 



dx 



dx 



24 



dx 



d^ 

de 



dx 



j_d^ 

2Ad^ 



dx 

"d^ 



Nesterenko recognized that there is a simple solution to this nonlinear problem: 



(40) 



(41) 



where Aq, m and a are constants. Substitution of this solution into the propagating equation leads to the following 
values: 



6(n- 2) 
n(n — 1) 



2\ 1/2 



CO = (2/71)1/2 



(42) 



An additional assumption is introduced at this point: a solitary wave is "constructed" by retaining this solution over 
one period, < a{k — cot) < tt, and setting dx/d^ equal to zero outside of this range. Note that this solution does 
not satisfy the velocity pulse initial condition, but rather it describes the solution that the system presumably settles 
into after a short initial transient, an assumption that is supported by numerical simulation results 0. For spherical 
granules 



dx 
"dk 



(fc - Cot). 



(43) 



If the initial velocity impulse is unity, then the initial total energy (all kinetic) is 1/2. The solitary wave H41|l . which 
describes both the potential and kinetic energy of the system once it settles, is assumed to contain essentially all of 
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this initial energy (the numerical simulation results confirm that the energy of the solitary wave is 99.7% of the initial 
energy 0). Because the potential is nonlinear, the potential and kinetic energies are no longer equal, but one can 
use the generalized equipartition theorem to calculate the average contribution of each. One finds that the ratio 
is K/U — ^, so that K — 2(n+2) ■ ^i^ice the velocity is 



the total kinetic energy in the pulse is 



where 



i(0 = co^osin— (44) 



^-4^^(r-T)-.T7r^. (45) 



2a 2(n + 2)' 



sin'M0 = 2' Y^fl)- 
The resulting pulse speed cq and pulse amplitude for ?i = 5/2 then are 

Co = 0.836, Ao = 0.765. (47) 

The numerical results of HSJ give cq — 0.84. 

Contrary to the n — 2 case, here there is almost no backscattering Except for the first two or three granules 
that are slightly scattered backwards, the granules in the pulse are simply displaced by a constant amount and come 
to rest once the pulse passes. The total backward momentum is thus extremely small and finite. Beyond the first 
two or three granules, the forward moving pulse here retains its shape and amplitude and is therefore essentially 
conservative with respect to both energy and momentum. This is to be contrasted with the fact that for cylindrical 
granules, every granule acquires a backward momentum. 



B. Granules With Friction - Theory 



When friction is included, Eq. H39|l is modified by the addition of a frictional term: 



at2 



dx d 
^~di ^ dk 



dk 



n-l 



n — 1 
24 



dx 
dk 



^ d^x 


1 d^ 


'(--)""' 


dk^ 


^ 24 9fc3 


\ dk) 



(48) 



For small values of 7, we expect the solution (|41|l to be modified in two ways. Firstly, as we did with the cylindrical 
granules, we must take into account the overall decay of the energy of the pulse, being mindful of the fact that kinetic 
and potential energies are not equal when n > 2. Secondly, since the speed of the pulse depends on its amplitude (and 
hence on its total energy), we must include the fact that the pulse speed decreases with time. We assume a solution 
of the form 



dx 



= yl(i)sin'«-2 a^{t), 



where 



and 



^(fc,t) = fc - / c{t)dt 
Jo 



n 



(49) 



(50) 



(51) 



Note that this form supposes that the width ir/a of the pulse is not changed by friction j2^ . 

The decay of A{t) [and hence of c{t)] can be determined by assuming that the pulse energy decays as e"^"'''* and 
choosing the constant u so that Eq. (|48|l is satisfied to first order in 7. If the pulse energy decays as e"^"'*'*, then the 
pulse velocity decays as e"""*"*. For the pulse velocity we have, aside from its ^ dependence, 



dx 
'dt 



c{t)A{t) ~ 



-u-yt 



(52) 
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from which it follows that 

A{t) = Aoe-^^\ c{t) = coe-^^^* (53) 

and therefore 

ak,t)^k-co " (l-e-^-^*). (54) 
uj(n — 2) \ / 

To determine u we note that Eq. H49|) implies that x{t,£^) — A{t)F{^), where the form of F is unimportant for the 
moment except that it depends only on ^ and not separately on t. Therefore 

Similarly we find 



dt dt \dtJ\d^J \A dt J \dt J 



d^x fld^A\ ^ f I dA\ f di\ f dx\ fd^^\fdx\ f d^V f d 



^2 



dx^ \Adt^ J'^ '^\Adt)\dt)\d^)^\dty\d^)^\dt) \de'' ^^^^ 

Substitution of Eqs. and (jS^ into (|^ and gives 

dx 2 2(n-2)u ..t d^x ^ ( , {n — 2)u 



W^^Tt^ + 7 ( - + '—r^ " 1 ) coe-^-* + 0(y ). (57) 



One must then choose 



(58) 



n + 2 

to force the 0(7) contribution to vanish. We are then left with 

For 7i < (n + 2)/2(n - 2), substitution into Eq. (gHl leads again to Eq.lgUJl to 0(7^). 

We thus conclude that to 0(7^) and for times 7i ^ (n + 2)/2{n — 2) (which means essentially the entire lifetime of 
the pulse, see next subsection) the solution for the chain of spherical granules subject to weak friction is as assumed 
in Eqs. and ISUJ) with 

Ait) = Aoe~^^\ c{t) = X -A7^e~^'". (60) 

V n 

The shape of the pulse is constant and the same as in the frictionless case. The width remains constant in time, the 
overall energy in the pulse decays as e~"^'^'^, and the pulse velocity as e~"T2^*. For n = 5/2 these decays go as 
g-(io/9)7t g-(5/9)7t respectively. It is an interesting sideline to note the increase in the friction-induced decay 
rate of the velocity or energy of the compression pulse with increasing n; for n — » 00 it is the same as that of a single 
particle traveling in a viscous medium. 

A remarkable difference between the frictionless and frictional chains of spherical granules lies in the backward 
scattering. In the frictionless case we noted that there is almost no backward scattering, and the very small amount 
of it appears only in the first three or so particles. With friction, each particle is scattered backward as the pulse 
leaves, more like the situation in the cylindrical granule case. Following the reasoning implemented in that case we 
first calculate the forward momentum of the pulse, 

P. E .(M). E (!)(§ 

x{k,t)>0 x{k,t)>0 ^ ^ 



^-^ a \n — 2 J 

4 1 (-^) e-*-^*, (61) 



x(k,t)>0 

(2/n)i/2 
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0.25 




1/n 



FIG. 5: Backscattering velocity at the moment of ejection as a function of the power n in the Hertz potential. The dashed line 
is the asymptote Eq. H66|l . which in this representation is independent of 7 and of t. The dark curve is the full result Eq. H()5|l 
for 7t = 1. 
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FIG. 6: Backscattering velocity at the moment of ejection as a function of the friction parameter for spherical granules at 
different times as follows: t = 10 (thick), 15 (dotted), 20 (dashed), 40 (dot-dashed), and 100 (thin). 




t 



FIG. 7: Total energy decay as a function of time. Within a time range of 0(107"^) the energy decays exponentially 
as e"'^"''^''^', drawn as lines. The symbols are numerical simulation results. From top to bottom the curves are for 
7 = 0.010, 0.012, 0.014, 0.016, 0.018, and 0.020. 
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where /(/) is given in Eq. (|46|l . The rate of change of the forward momentum then is 



{n + 2) " \n-2^ 

On the other hand, the loss of momentum of the last particle in the pulse as the pulse moves across it is [cf. Eq. (|32|l ] 



a \n — 2 



a \n — 2 

The momentum transferred to the last particle as it is ejected from the pulse is therefore 



v,^-JL + s = -2j /° j f^)e-^^* (64) 



c{t) 'a{n + 2) \n-2 

(in Eqs. p5|) and (|33|l the factor c{t) does not appear explicitly because it is equal to unity). For spherical granules 

5 ttAo „4^j 



'"' = ~^\J2~^^ - ^^^^ 

Every particle in the chain is therefore ejected with a (exponentially decreasing) backward momentum, in contrast 
to the frictionless chain. The loss of energy in the pulse due to friction is not balanced by a sufficiently large loss of 
momentum. This relative increase of pulse momentum must be balanced by the momentum carried by the ejected 
granules. 

From Eq. I|64|l we can deduce two additional interesting results. One is the dependence of the backward velocity 
on the power n of the potential of interaction. We find that Vh decreases monotonically with n, and for large n, 

n \ J 

independently of time to leading order. The approach to this asymptotic result is shown in Fig. [3 where we present 
ff,/7 as a function of 1/n. The other is the nonmonotonic dependence of Uf, on the damping, a result already found 
in our two-granule analysis j21.] and illustrated in Fig. |^ For spherical granules we show this dependence of Vb on 
damping at different times. This qualitative behavior persists for other values oi n > 2: as n increases the position of 
the minimum is shifted to larger 7 and the absolute value of the minimum decreases. 

Equation (|49|l with (|60l) and Eq. H64() are the principal results of this subsection. They establish analytic expressions 
for the pulse and for the backscattered momentum that will be checked against numerical simulations in the following 
subsection. 



C. Numerical Simulations 



First we note that even in the frictionless case the pulse velocity is lower here (cq = 0.836) than in the cylindrical 
case (unit velocity), a result obtained in the numerical simulations of HSJ as well as in our theory, cf. Eq. H47() . With 
friction, the pulse speed decreases as the pulse loses energy. 

As with cylindrical granules, for times t < 107^^ the pulse has the same shape in the presence of friction as in the 
frictionless case except for the overall exponential decay, as illustrated in Fig. [7| Note that the decay is more rapid 
than in the case of cylindrical granules. 

Contrary to the cylindrical granule case, the pulse also slows down as its energy decreases. This behavior continues 
until the energy begins to decrease abruptly (more rapidly than exponentially) at a time ~ 8.37"^. At a time 
~ 16.27"'^ the backscattered energy becomes greater than the pulse energy and at the same time the pulse stops 
moving. In Fig.|Slwe show kmax, the granules with the maximum velocity in the pulse. The symbols are the results 
obtained from numerical simulations, and the lines are [cf. Eq. H54|) ] 
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FIG. 8: Grain with the maximum velocity as a function of time. The symbols are simulation results, and the lines are Eq. (I67I I. 
From top to bottom 7 = 0, 0.001, 0.005, 0.01, 0.012, 0.014, 0.016, 0.018, 0.02, 0.03. The 7 = curve is k^ax = (7r/2a) + cot. 
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FIG. 9: Scaled velocity pulse as a function of scaled position for different friction coefficients. 



where the last expression holds expUcitly for n = 5/2. The value used for cq is as given in Eq. (|47|) . The agreement 
is clearly excellent. 

To check our prediction for the effect of friction on the propagating pulse we present the scaled velocity pulse as a 
function of scaled position for different friction coefficients in Fig. |^ The agreement here is less satisfactory than in 
the cylindrical granule case. One difficulty is the fact that the analytic theory is a continuum approximation while 
the chain is discrete, and here only a very small number of granules are actually moving forward at any one time, 
that is, the pulse is very narrow. A second difficulty is that in addition to the smooth envelope that the continuum 
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FIG. 10: Backscattering velocity at the moment of ejection as a function of time for different friction coefficients. The 
symbols are simulation results and the lines are Eq. (|65^ . From top to bottom on the right side of the figure 7 = 
0.010, 0.012, 0.014, 0.016, 0.018, 0.020. 



theory attempts to capture, the narrow pulse actually experiences small amplitude oscillations as it moves forward. 
Our reported values in Fig. include values that might fall anywhere within these oscillations. It is nevertheless clear 
that the theory captures the qualitative features of the pulse. In particular, we point to the increasing backscattering 
with increasing friction that can be seen in these figures. This is a qualitative difference between the frictionless chain 
(where only about three granules backscatter slightly) to the chain with friction, where some backscattering occurs at 
each granule as the pulse passes by. In Fig. llOl the symbols are the simulation results and the lines represent Eq. Ht)5|l . 
The agreement is clearly excellent. 



IV. CONCLUSIONS 



In this work we have studied the dynamics of an initial velocity impulse in a chain of granules that interact only 
when in contact, that is, they experience only a repulsive potential. Our interest has been in establishing the effects 
of hydrodynamic friction on these dynamics. 

First we analyzed a chain of cylindrical granules, that is, with a power law repulsive potential with exponent n = 2 
(half a harmonic potential). We presented the frictionless case, and organized existing information in a particular way 
to clarify the effects of discreteness and of the absence of restoring forces on these results. We were also able to obtain 
a number of results analytically that had previously only been obtained numerically 1^. In this chain the impulse 
travels at unit velocity as a spreading pulse (as t^/"^) in which the maximum displacement progressively grows with 
time (as t^^^). While this traveling pulse carries most of the initial energy, conservation of momentum requires that 
there be backscattering of each granule as the pulse passes by. The chain thus continually undergoes fragmentation. 

We then generalized these results in the presence of friction and found that the principal effect of weak friction is an 
overall exponential decay of the energy. The pulse still moves at unit velocity, still spreads as t^^^, and the maximum 
displacement now varies as t^/'^i"'''*/^ throughout its lifetime. There is again backscattering of each granule as the 
pulse passes by. An interesting and unexpected effect of friction is that the velocity of the backscattered particles at 
the moment of ejection is greater than the velocity in the frictionless chain '2l]. The backscattered particles also slow 
down due to friction, but this chain, too, undergoes continual fragmentation. We supported our results via numerical 
simulations. 

Next we analyzed a chain of granules with a power law repulsive potential with exponent n > 2, with special 
attention to spherical granules (n = 5/2). We reviewed Nesterenko's theory for the frictionless case, which leads to 
an essentially conservative pulse of constant width determined by the power n, traveling down the chain at a velocity 
that depends on the energy of the pulse. This velocity is lower than that of the spreading pulse in the cylindrical 
granule case. Here again we obtained some results analytically that had previously been reported numerically |^. 
Contrary to the n = 2 case there is essentially no backscattering (fragmentation) in this system: only the first two or 
three granules acquire a very small backward momentum as the pulse passes over. 

The generalization of these results in the presence of friction are more complicated because an overall decay of the 
energy causes the pulse to slow down as it moves. We found that the solution is one in which the overall shape of the 
pulse as well as its width remain unchanged, the energy decay is exponential, as is the decrease in the displacement 
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pulse amplitude and the pulse velocity. The exponential decay factors are fixed by the power n of the potential. Most 
dramatically, we found that in the presence of friction there is now backscattering of each granule as the pulse passes 
by, so that this chain experiences fragmentation. The velocity of the backscattered granules of course also decreases 
exponentially. We supported these results for the case of spherical granules via numerical simulations. 

A number of interesting problems immediately come to mind as a possible extension of this work. Among them is 
consideration of the effect of mass disorder and/or frictional disorder in the chains and of mass tapering. Work along 
these directions is in progress |3C| . 
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